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This  final  report  summarizes  accomplishments  in  two  areas  of  uncer¬ 
tainty  quantification  and  computational  probability. 

1  Methods  for  modeling  and  integrating  different  types  of 
uncertainty 

1.1  Problem  formulation  and  issues 

Uncertainty  quantification  refers  to  a  broad  set  of  techniques  for  understand¬ 
ing  the  impact  of  uncertainties  in  complicated  mechanical  and  physical  sys¬ 
tems.  In  this  context  ‘‘uncertainty”  can  take  on  many  meanings.  Aleatoric 
uncertainty  refers  to  inherent  uncertainty  due  to  stochastic  or  probabilistic 
variability.  This  type  of  uncertainty  is  irreducible  in  that  there  will  always 
be  positive  variance  since  the  underlying  variables  are  truly  random.  Epis- 
temic  uncertainty  refers  to  limited  knowledge  we  may  have  about  the  model 
or  system.  This  type  of  uncertainty  is  reducible  in  that  if  we  have  more  infor¬ 
mation,  e.g.,  take  more  measurements,  then  this  type  of  uncertainty  can  be 
reduced.  For  many  problems  where  uncertainty  quantification  is  important, 
the  acquisition  of  data  is  difficult  or  expensive.  The  epistemic  uncertainty 
cannot  be  removed  entirely,  and  so  one  needs  modeling  and  computational 
techniques  which  can  also  accommodate  this  form  of  uncertainty. 

However,  in  most  applications  it  is  (perhaps  implicitly)  assumed  that 
epistemic  uncertainty  can  be  modeled  by  aleatoric  uncertainty.  One  reason 
is  purely  based  on  convenience  in  that,  at  least  as  they  have  been  developed 
to  date,  most  computational  techniques  (e.g.,  polynomial  chaos  and  Monte 
Carlo)  are  based  on  the  assumption  that  the  user  can  identify  some  “ap¬ 
propriate”  distribution  for  each  uncertain  aspect  of  the  system,  regardless 
of  the  type  of  uncertainty,  aleatoric  or  epistemic.  If  one  is  interested  in  just 
basic  qualitative  properties  of  the  system  then  this  may  not  be  a  central 
issue,  since  virtually  any  model  of  uncertainty  will  give  information  on  the 
sensitivities  of  the  system.  However,  when  the  intended  use  of  uncertainty 
quantification  is  for  regulatory  assessment  or  some  other  application  where 
performance  measures  are  sensitive  to  distributional  assumptions,  the  issue 
becomes  much  more  important,  and  one  should  carefully  distinguish  how 
one  accounts  for  the  two  types  of  uncertainty. 

Hence  important  issues  are:  (i)  how  should  one  model  uncertainties  that 
are  not  of  the  aleatoric  type,  and  (ii)  can  one  work  with  the  resulting  for¬ 
mulation  from  a  numerical  perspective. 
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1.2  Results  and  key  findings 


In  the  paper  [1]  we  develop  an  approach  that  (i)  logically  distinguishes  those 
aspects  of  uncertainty  that  are  treated  as  stochastic  variability  from  other 
forms  of  uncertainty,  (ii)  in  cases  where  a  stochastic  model  is  theoretically 
valid  but  for  which  determination  of  the  distribution  is  not  practical,  gives 
bounds  for  performance  measures  that  are  valid  for  explicitly  identified  fam¬ 
ilies  of  distributions,  and  (iii)  is  computationally  feasible  if  ordinary  uncer¬ 
tainty  propagation  is  feasible.  The  different  forms  of  uncertainty  that  are 
covered  by  the  formulation  include:  (a)  aleatoric  with  known  distribution; 
(b)  aleatoric  with  partly  known  distribution  (mingled  aleatoric  and  epis- 
temic);  (c)  epistemic  for  which  one  is  willing  to  model  by  a  family  of  aleatoric 
uncertainties,  and  (d)  epistemic  where  one  is  only  willing  to  place  bounds 
on  the  uncertainties. 

Suppose  that  random  variables  with  a  known  distribution  will  take  val¬ 
ues  in  a  Polish  space  X.  Variables  whose  distribution  is  not  known  or  are 
otherwise  of  the  epistemic  variety  take  values  in  the  space  y.  Let  the  per¬ 
formance  measure  of  interest  for  some  given  problem  is  assumed  to  be  of  the 
form 


F(x,y)'y(dy)^(dx), 


where  (i  (resp.,  7)  is  a  probability  measure  on  X  (resp.,  y).  If  X  and  Y 
are  independent  random  variables  with  distributions  ^  and  7,  then  F(V,  Y) 
represents  both  the  performance  measure  (e.g.,  a  second  moment)  as  well 
as  the  underlying  physical  or  mechanical  system  that  maps  these  aleatoric 
and  epistemic  inputs  into  outputs. 

It  is  known  that  risk-sensitive  performance  measures  can  be  used  to 
produce  performance  bounds  that  are  robust  with  respect  to  the  underlying 
distributions.  The  standard  risk-sensitive  performance  measure  is 


Ac  =  —  log  f  f  ecF(x’y')') (dy)y(dx).- 

c  Jx  Jy 

Neither  of  these  measures  differentiate  the  variables  according  to  type  (aleatoric 
or  epistemic)  and  given  that  the  performance  measure  of  interest  is  actually 
F,  the  use  of  a  risk-sensitive  version  of  the  cost  is  not  well  motivated  for 
the  aleatoric  variables.  Indeed,  use  of  this  measure  will  give  bounds  that 
are  robust  with  respect  to  variations  on  a  distribution  that  is  known,  and 
obviously  such  bounds  will  not  be  as  tight  as  possible. 

In  [1]  two  hybrid  risk-sensitive  measures  are  introduced,  as  well  as  vari¬ 
ations.  To  simplify  we  mention  only  the  theory  developed  for  the  (more 
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useful)  form 


Aj  =  -  log  [ 

C  Jy 

Using  the  relative  entropy  representation  for  exponential  integrals,  it  follows 
that  for  any  distribution  0(dy) 

[  [  F(x,  y)[*(dx)0(dy)  <  -R  ( 6{dy )  ||'y(rfy) )  +  A*. 

JyJx  c 

This  gives  a  bound  on  the  performance  measure  for  an  arbitrary  distribution 
on  Y,  but  with  the  distribution  on  X  equal  to  the  known  true  distribution. 
The  distributions  thus  play  very  different  roles.  In  particular,  we  think  of  7 
as  a  nominal  distribution  of  K,  which  should  be  distinguished  from  a  pos¬ 
sible  true  distribution.  The  risk  sensitive  functional  A*,  whose  numerical 
evaluation  can  be  carried  out  by  a  variety  of  methods  (in  [1]  we  use  poly¬ 
nomial  chaos),  is  based  on  the  nominal  distribution.  Through  the  relative 
entropy  duality,  it  yield  various  bounds  (depending  on  c)  on  a  families  of 
distributions,  with  the  relative  entropy  distance  the  key  metric. 

Suppose  that  a  bound  on  performance  over  a  specific  family  of  distri¬ 
butions  is  needed.  Let  R*  denote  the  maximum  of  relative  entropy  with 
respect  to  the  nominal  model  over  this  family  that  we  wish  to  allow.  Then 
the  tightest  possible  bound  is  obtained  by  minimizing 

A  *  +  -R* 

C 

over  c  >  0.  Thus  given  a  prescribed  uncertainty  in  the  epistemic  variables, 
one  can  compute  a  bound  on  the  performance  over  all  distributions  in  the 
family  and  which  holds  with  equality  for  at  least  one.  This  optimizing  c, 
which  exists  and  is  unique,  can  be  computed  using  the  same  techniques  used 
to  compute  the  performance  measure  itself.  We  show  in  [1]  that  this  function 
has  only  one  local  minimum  over  c  £  (0,oo],  and  thus  the  global  minimum 
is  easy  to  compute. 

In  [1]  examples  of  various  types  and  numerical  data  is  presented,  as 
well  as  the  corresponding  theory  where  all  that  is  assumed  regarding  the 
epistemic  variables  is  that  they  are  constrained  to  some  given  set.  The  use  of 
polynomial  chaos  methods  for  the  evaluation  of  the  risk-sensitive  functionals 
is  developed,  and  explicit  relative  entropy  distances  for  common  families  of 
distributions  are  listed. 

Two  additional  papers  on  this  topic  are  in  preparation,  both  with  Kenny 
Chowdhary,  a  graduate  student  in  the  department.  One  is  concerned  with 
optimization  in  the  presence  of  both  epistemic  and  aleatoric  uncertainties, 
and  the  second  is  concerned  with  estimation  in  this  same  framework. 
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2  Computational  methods  for  rare  event  problems 

2.1  Problem  formulation  and  issues 

There  is  significant  interest  in  the  application  of  uncertainty  quantification  to 
problems  with  small  probabilities  (rare  events).  For  example,  it  is  a  research 
focus  for  the  SAMSI/Sandia  Summer  School  on  Uncertainty  Quantification 
in  2011.  Though  these  probabilities  may  be  small,  they  are  often  critical 
measures  of  system  performance  and  one  needs  reasonably  reliable  numerical 
methods.  Unfortunately,  standard  numerical  schemes  are  not  at  all  reliable 
for  problems  with  rare  events. 

This  part  of  the  research  project  was  concerned  with  developing  efficient 
Monte  Carlo  algorithms  for  rare  event  simulation  and  the  associated  large 
deviations  theory.  There  are  two  methods  currently  in  use.  One  is  based 
on  simulating  according  to  a  different  distribution  and  then  correcting  for 
any  induced  bias  via  the  likelihood  ratio  (importance  sampling).  The  key 
question  here  is  how  to  select  the  new  sampling  distribution.  The  second 
method  simulates  a  branching  process,  i.e.,  collection  of  particles  that  can 
split  according  to  certain  rules  to  form  new  particles,  each  of  which  behaves 
like  the  original  particle  or  process.  The  splitting  rules  are  designed  to 
make  the  rare  event  likely  for  at  least  one  of  the  descendent  particles,  and 
the  estimator  is  the  ratio  of  number  of  particles  for  which  the  rare  outcome 
is  observed  to  the  total  number  of  descendents.  The  key  question  here 
is  what  should  trigger  a  split  and,  given  that  a  split  occurs,  the  number 
of  descendents.  Most  of  the  literature  on  these  methods  features  schemes 
based  on  heuristic  design,  with  little  or  no  analysis.  However,  the  design 
problem  with  both  methods  is  subtle,  and  reasonable  looking  schemes  can 
perform  quite  poorly.  Indeed,  simulations  based  on  improperly  designed 
schemes  could  be  highly  misleading. 

A  second  class  of  problems  considers  the  numerical  approximation  of  the 
invariant  distribution  for  stochastic  systems  with  multiple  metastable  states 
using  the  occupation  measure  of  a  related  Markov  process.  Moving  from 
one  metastable  state  to  another  is  a  rare  event,  and  its  treatment  is  the  key 
question  in  the  design  of  efficient  Monte  Carlo  schemes.  There  are  many 
ad  hoc  algorithms  available.  However,  these  algorithms  do  not  always  work 
well  and  have  to  be  applied  with  some  care. 

2.2  Results  and  key  findings 

In  prior  work  we  demonstrated  that  the  design  of  a  reliable  and  high  per¬ 
forming  important  sampling  scheme  would  follow  if  one  could  construct  a 
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subsolution  to  a  related  Hamilton- Jacobi-Bellman  equation  (a  partial  dif¬ 
ferential  equation).  Moreover,  we  showed  that  the  existence  of  such  a  sub¬ 
solution  is  in  some  sense  necessary.  The  paper  '4]  considers  the  analysis 
of  the  “weighted-serve-the-longer-queue”  policy.  Such  service  policies  are 
common,  and  in  particular  a  variant  is  frequently  used  in  wireless  com¬ 
munication.  The  goal  of  the  importance  sampling  scheme  is  to  accurately 
estimate  probabilities  associated  with  buffer  overflows  and  delays,  which  are 
critical  performance  measures.  The  model  is  fairly  complex,  since  the  policy 
introduces  a  kind  of  discontinuous  behavior  into  the  stochastic  evolution. 
We  showed  how  the  approach  based  on  subsolutions  to  a  related  Hamilton- 
Jac obi- Bellman  equation  developed  in  previous  papers  could  be  adapted  to 
deal  with  such  complicated  process  dynamics. 

The  paper  [2]  analyzes  a  number  of  branching  type  schemes,  including 
RESTART  (REpetitive  Simulation  Trials  After  Reaching  Thresholds)  and 
DPR  (Direct  Probability  Redistribution).  It  is  established  that  the  design  of 
a  stable  (namely,  controlled  number  of  particles)  and  asymptotically  optimal 
(namely,  tightly  controlled  variance)  splitting  algorithms  can  be  achieved  by 
constructing  suitable  viscosity  subsolutions  to  a  Hamilton- Jacobi-Bellman 
(HJB)  equation.  This  HJB  equation  is  in  fact  the  same  one  in  importance 
sampling  analysis.  This  is  a  useful  theoretical  result,  since  it  indicates  that 
the  construction  of  subsolutions  is  in  any  case  a  central  aspect  to  solving 
the  numerical  problem. 

A  much  more  complex  algorithm  can  be  based  on  what  is  called  an  in¬ 
teracting  particle  system .  The  paper  [3]  provides  the  first  rigorous  analysis 
of  the  performance  of  this  class  of  algorithms  in  the  small  probability  limit. 
Owing  to  the  complexities  of  the  algorithm,  the  analysis  is  limited  to  di¬ 
mension  one.  However,  the  results  support  some  of  the  claims  that  have 
been  made  concerning  these  algorithms  (but  based  only  on  numerical  evi¬ 
dence),  and  in  particular  that,  at  least  for  one  dimensional  problems,  they 
are  less  sensitive  to  the  details  of  the  underlying  distributions  than  compet¬ 
ing  schemes. 

The  paper  [5]  considers  the  problem  of  approximating  stationary  distri¬ 
butions  of  a  Markov  chain  by  simulation.  Our  initial  goal  was  to  use  large 
deviation  ideas  to  choose  design  parameters  in  an  existing  scheme  known  as 
parallel  tempering.  Parallel  tempering  (also  known  as  replica  exchange  sam¬ 
pling)  is  a  standard  method  for  simulating  complex  systems,  and  is  used  in 
many  commercial  software  packages.  In  this  algorithm  simulations  are  con¬ 
ducted  in  parallel  for  a  family  of  Markov  chains  indexed  by  a  “temperature” 
parameter,  and  the  key  improvement  over  standard  Monte  Carlo  is  a  swap 
mechanism  that  exchanges  configurations  between  these  parallel  simulations 
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at  a  given  rate.  The  mechanism  is  designed  to  allow  the  low  temperature 
system  to  escape  from  deep  local  energy  minima  where  it  might  otherwise  be 
trapped,  via  those  swaps  with  the  higher  temperature  components.  Based 
on  large  deviation  theory,  we  have  argued  that  the  rate  of  convergence  of 
the  empirical  measure  is  a  monotone  increasing  function  of  the  swap  rate. 
This  suggests  that  one  should  raise  the  swap  frequency  in  order  to  improve 
efficiency,  but  this  is  eventually  counter-productive  since  eventually  most 
of  the  computational  effort  is  directed  towards  swapping  and  little  towards 
moving  the  process  dynamics.  However,  it  turns  out  that  one  can  construct 
a  simulation  scheme  that  is  equivalent  to  the  limit  of  the  parallel  tempering 
schemes  in  the  sense  of  distributions,  but  which  involves  no  swapping  at 
all.  With  this  scheme,  wdiich  w'e  call  infinite  swapping  (INS),  the  effect  of 
the  swapping  is  captured  by  a  collection  of  weights  that  influence  both  the 
dynamics  and  the  empirical  measure. 

While  the  infinite  swapping  scheme  optimizes  the  convergence  rate  as 
described  above,  it  has  practical  limitations  in  implementation  due  to  the 
complexity  of  the  weights  when  the  number  of  temperatures  is  large  (>7). 
Complex  problems  in  often  involve  scores  of  temperatures,  and  so  it  was 
critical  to  overcome  this  limitation.  We  have  developed  an  approximation 
to  the  full  infinite  swapping  which  is  based  on  alternating  between  partial 
infinite  swapping  (PINS)  algorithms,  wffiich  can  be  shown  to  approximate 
(theoretically  and  practically)  the  INS  scheme.  The  mathematical  theory  for 
the  INS  and  PINS  is  developed  in  [5].  Numerical  studies  on  fairly  complex 
Lennar d- Jones  systems  (very  challenging  benchmark  problems  from  chem¬ 
istry)  have  been  conducted.  Improvements  of  three  orders  of  magnitude 
in  performance  over  conventional  parallel  tempering  were  observed  at  an 
increased  computational  cost  of  5-15%. 
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